﻿#pragma once

#include "../../Utility/Util/utility.hpp"
#include "../../Utility/Log/Log.hpp"
#include "../../Utility/Static/Mathf.h"
#include "../../Utility/Static/Matrixf.h"
#include "../../Utility/Data/Position/CRDFactory.hpp"
#include "../../Utility/Data/Time/TimeFactory.hpp"
#include "../../Utility/Data/NavMessageData/GNSSMesData.hpp"

namespace ns_Calculate
{
    using namespace ns_Time;
    using namespace ns_Position;
    using namespace ns_GNSSMesData;
    using namespace ns_Util;
    using namespace ns_Log;

    class IonoCorrect
    {
    public:
        IonoCorrect()
        {
        }

        // 获得Klobuchar模型的电离层改正参数，单位为秒
        static double GetKlobuchar(GNSSType type, const CRDFactory::CRDPtr &rr, const CRDFactory::CRDPtr &rs,
                                   const GNSSMesHeader::Ptr &header, const TimeFactory::TimePtr &curTime)
        {
            CRDFactory::CRDPtr _star = CRDFactory::MakeCRDPtr(CRDFactory::ValueType::TopoCentricPolarCRDSys, CartesianCRD());
            CRDFactory::SwitchTo(CRDFactory::ValueType::TopoCentricPolarCRDSys, rs, &_star, rr);
            auto star = std::dynamic_pointer_cast<TopoCentricPolarCRD>(_star);

            CRDFactory::CRDPtr _geoRr = CRDFactory::MakeCRDPtr(CRDFactory::ValueType::GeodeticCRDSys);
            CRDFactory::SwitchTo(CRDFactory::ValueType::GeodeticCRDSys, rr, &_geoRr, nullptr);
            auto geoRr = std::dynamic_pointer_cast<GeodeticCRD>(_geoRr);

            TimeFactory::TimePtr _curTimeOfGPSTime = TimeFactory::MakeTimePtr(TimeFactory::ValueType::GPSTimeSys);
            TimeFactory::SwitchTo(TimeFactory::ValueType::GPSTimeSys, curTime, &_curTimeOfGPSTime);
            auto curTimeOfGPSTime = std::dynamic_pointer_cast<GPSTime>(_curTimeOfGPSTime);

            switch (type)
            {
            case GNSSType::GPS:
                return GPSKlobucharCalculate(*star, geoRr.get(), *header, *curTimeOfGPSTime);

            case GNSSType::BDS:
                return BDSKlobucharCalculate(*star, geoRr.get(), *header, *curTimeOfGPSTime);

            default:
                Log(Debug) << "No such Correct Model" << '\n';
            }

            return -1;
        }

    private:
        static double GPSKlobucharCalculate(const TopoCentricPolarCRD &input, GeodeticCRD *geoRr, const GNSSMesHeader &header, const GPSTime &curTime)
        {
            double centralAngle = 0.0137 / ((input.GetElevation() / Mathf::Pi) + 0.11) - 0.022;

            // 电离层穿刺点的纬度
            double FAIi = geoRr->GetB() / 180 + centralAngle * Mathf::Cos(input.GetAzimuth());
            if (FAIi > 0.416)
                FAIi = 0.416;
            if (FAIi < -0.416)
                FAIi = -0.416;

            // 电离层穿刺点经度
            double Lambdai = geoRr->GetL() / 180 + (centralAngle * Mathf::Sin(input.GetAzimuth())) / Mathf::Cos(FAIi * Mathf::Pi);

            // 电离层穿刺点的地磁纬度
            double FAIm = FAIi + 0.064 * Mathf::Cos((Lambdai - 1.617) * Mathf::Pi);

            // GPS周内秒
            double timeTmp = 43200 * Lambdai + curTime.GetSecond();
            double t = timeTmp - Mathf::Floor(timeTmp / 86400) * 86400;

            // 电离层延迟的幅度
            double Ai = header.ionAlpha[0] + FAIm * (header.ionAlpha[1] + FAIm * (header.ionAlpha[2] + FAIm * header.ionAlpha[3]));
            if (Ai < 0)
                Ai = 0;

            // 电离层延迟的周期
            double cyclei = header.ionBeta[0] + FAIm * (header.ionBeta[1] + FAIm * (header.ionBeta[2] + FAIm * header.ionBeta[3]));
            if (cyclei < 72000)
                cyclei = 72000;

            // 电离层延迟的相位
            double Xi = 2 * Mathf::Pi * (t - 50400) / cyclei;

            // 倾斜因子
            double F = 1.0 + 16.0 * Mathf::Pow(0.53 - input.GetElevation() / Mathf::Pi, 3);

            // 计算结果
            if (Mathf::Abs(Xi) >= 1.57)
            {
                return 5e-9 * F * Mathf::c;
            }
            else
            {
                return (5E-9 + Ai * (1 - 0.5 * Xi * Xi + Mathf::Pow(Xi, 4) / 24.0)) * F * Mathf::c;
            }
        }

        static double BDSKlobucharCalculate(const TopoCentricPolarCRD &input, GeodeticCRD *geoRr, const GNSSMesHeader &header, const GPSTime &curTime)
        {
            double centralAngle = 0.0137 / ((input.GetElevation() / Mathf::Pi) + 0.11) - 0.022;

            double R = 6378000; // 地球半径
            double h = 375000;  // 电离层单层高度

            // 用户地理经纬度
            double Faiu = geoRr->GetB() / 180;
            double Lambdau = geoRr->GetL() / 180;

            // 用户和穿刺点的地心张角
            double Fai = Mathf::Pi / 2 - input.GetElevation() - Mathf::Arcsin((R) / (R + h) * Mathf::Cos(input.GetElevation()));

            // 电离层穿刺点的地理经纬度
            double Faim = Mathf::Arcsin(Mathf::Sin(Faiu) * Mathf::Cos(Fai) + Mathf::Cos(Faiu) * Mathf::Sin(Fai) * Mathf::Cos(input.GetAzimuth()));
            double Lambdam = Lambdau + Mathf::Arcsin(Mathf::Sin(Fai) * Mathf::Sin(input.GetAzimuth()) / Mathf::Cos(Faim));

            // GPS周内秒
            double timeTmp = 43200 * Lambdam / Mathf::Pi + curTime.GetSecond();
            double t = timeTmp - Mathf::Floor(timeTmp / 86400) * 86400;

            // 白天电离层延迟余弦曲线的幅度
            double A2 = header.ionAlpha[0] + Mathf::Abs(Faim / Mathf::Pi) * (header.ionAlpha[1] + Mathf::Abs(Faim / Mathf::Pi) * (header.ionAlpha[2] + Mathf::Abs(Faim / Mathf::Pi) * header.ionAlpha[3]));
            if (A2 < 0)
                A2 = 0;

            // 余弦曲线的周期
            double A4 = header.ionBeta[0] + Mathf::Abs(Faim / Mathf::Pi) * (header.ionBeta[1] + Mathf::Abs(Faim / Mathf::Pi) * (header.ionBeta[2] + Mathf::Abs(Faim / Mathf::Pi) * header.ionBeta[3]));
            if (A4 > 172800)
                A4 = 172800;
            else if (A4 < 72000)
                A4 = 72000;

            // B1I信号的电离层垂直延迟改正
            double Izt;
            if (Mathf::Abs(t - 50400) < A4 / 4)
            {
                Izt = 5E-9 + A2 * Mathf::Cos(2 * Mathf::Pi * (t - 50400) / A4);
            }
            else
            {
                Izt = 5E-9;
            }

            // B1I信号传播路径上的电离层延迟
            double IB1It = 1 / Mathf::Sqrt(1 - Mathf::Pow(R / (R + h) * Mathf::Cos(input.GetElevation()), 2)) * Izt;

            return IB1It * Mathf::c;
        }
    };
}